if (rstudioapi::isAvailable()) {
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
}library(tidyverse)## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
## ✔ ggplot2 3.3.6 ✔ purrr 0.3.4
## ✔ tibble 3.1.8 ✔ dplyr 1.0.10
## ✔ tidyr 1.2.0 ✔ stringr 1.4.1
## ✔ readr 2.1.2 ✔ forcats 0.5.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
library(conflicted)
conflict_prefer("filter", "dplyr")## [conflicted] Will prefer dplyr::filter over any other package
conflict_prefer("select", "dplyr")## [conflicted] Will prefer dplyr::select over any other package
conflict_prefer("slice", "dplyr")## [conflicted] Will prefer dplyr::slice over any other package
conflict_prefer("rename", "dplyr")## [conflicted] Will prefer dplyr::rename over any other package
conflict_prefer('intersect', 'dplyr')## [conflicted] Will prefer dplyr::intersect over any other package
table with all measured Cq values
cq = read_tsv("../data/details/selected_circ_cq.txt")## Rows: 1560 Columns: 25
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (10): chr, strand, circ_id, circ_id_strand, tool, cell_line, FWD_primer,...
## dbl (15): start, end, FWD_pos, FWD_length, REV_pos, REV_length, FWD_TM, REV_...
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
cqcq = cq %>%
# qPCR_val is 'pass' for all circ that have at least one replicate > 10
mutate(qPCR_val = ifelse((Cq_min_untreated < 10) & (Cq_max_untreated < 10), 'fail', 'pass'),
qPCR_val = ifelse((is.na(Cq_min_untreated)) & (is.na(Cq_max_untreated)), 'fail', qPCR_val),
# qPCR_val_detail describes the reason why a circ fails this validation step if it does
qPCR_val_detail = ifelse((Cq_min_untreated < 10) & (Cq_max_untreated < 10), 'all_Cq_under_10', 'pass'),
qPCR_val_detail = ifelse(is.na(Cq_min_untreated) & is.na(Cq_max_untreated), 'all_NA', qPCR_val_detail))
cqstatistics
cq %>% count(qPCR_val)cq %>% count(qPCR_val_detail)cq = cq %>%
# calculate delta cq
# first make temp 47 for NAs
mutate(Cq_min_treated_tmp = Cq_min_treated,
Cq_min_treated = ifelse(is.na(Cq_min_treated), 47, Cq_min_treated),
Cq_diff = Cq_min_treated - Cq_max_untreated) %>%
mutate(# fail if diff > 3
RR_val = ifelse(Cq_diff > 3, "fail", 'pass'),
RR_val_detail = ifelse(RR_val == 'fail', "FP (Cq_diff > 3)", 'pass'),
# when Cq untreated is high => put NA
RR_val = ifelse(Cq_min_untreated > 32, NA, RR_val),
RR_val_detail = ifelse(Cq_min_untreated > 32, 'out_of_range', RR_val_detail),
# if qPCR already failed => RR also becomes 'fail'
RR_val = ifelse(qPCR_val == "fail", 'fail', RR_val),
RR_val_detail = ifelse(qPCR_val == "fail", "fail_qPCR", RR_val_detail)) %>%
# clean up NAs for 47
mutate(Cq_diff = ifelse(Cq_min_treated == 47, NA, Cq_diff),
Cq_min_treated = Cq_min_treated_tmp) %>%
select(-Cq_min_treated_tmp) %>%
#' to clean up dataframe: also remove Cq_diff for 'all_Cq_under_10' group and out_of_range group
mutate(Cq_diff = ifelse(qPCR_val_detail == 'all_Cq_under_10', NA, Cq_diff),
Cq_diff = ifelse(RR_val_detail == 'out_of_range', NA, Cq_diff))
cqstatistics
cq %>% count(RR_val)cq %>% count(RR_val_detail)cq_amp = read_tsv('../data/details/amp_seq.txt')## Rows: 1216 Columns: 5
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (4): circ_id, cell_line, amp_seq_val, amp_seq_val_detail
## dbl (1): perc_on_target
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
statistics
cq_amp %>% count(amp_seq_val)cq_amp %>% count(amp_seq_val_detail)add to cq dataframe
cq = cq %>%
left_join(cq_amp %>% unique()) %>%
mutate(amp_seq_val_detail = ifelse(is.na(amp_seq_val_detail), "not_included", amp_seq_val_detail))## Joining, by = c("circ_id", "cell_line")
statistics
cq %>% count(amp_seq_val)cq %>% count(amp_seq_val_detail)cq = cq %>%
mutate(RR_val_tmp = ifelse(RR_val_detail == 'out_of_range', 'fail', RR_val)) %>% # when circ cannot be measured for RR => count as fail
mutate(compound_val = ifelse(qPCR_val == "pass" & RR_val_tmp == "pass" & amp_seq_val == "pass", "pass", "NA"),
compound_val = ifelse(qPCR_val == 'fail' | RR_val_tmp == "fail" | amp_seq_val == 'fail', 'fail', compound_val),
compound_val = ifelse(is.na(amp_seq_val), NA, compound_val)) %>%
select(-RR_val_tmp)statistics
cq %>% count(compound_val)all_circ = read_tsv("../data/Supplementary_Table_2_all_circRNAs.txt")## Rows: 1137099 Columns: 23
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (16): chr, strand, cell_line, tool, circ_id, circ_id_strand, count_group...
## dbl (7): start, end, BSJ_count, n_detected, n_db, estim_len_in, BSJ_count_m...
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
all_circcq = cq %>%
left_join(all_circ) ## Joining, by = c("chr", "start", "end", "strand", "circ_id", "circ_id_strand",
## "tool", "cell_line", "BSJ_count", "count_group")
cq %>%
write_tsv('../data/Supplementary_Table_3_selected_circRNAs.txt')perc_val = cq %>%
select(tool, circ_id, cell_line, count_group, qPCR_val, RR_val, RR_val_detail, amp_seq_val, amp_seq_val_detail, compound_val) %>%
group_by(tool, count_group) %>%
summarise(nr_qPCR_total = n(),
nr_qPCR_fail = sum(qPCR_val == 'fail'),
nr_qPCR_val = sum(qPCR_val == 'pass'),
nr_RR_total = n() - sum(is.na(RR_val)), # here NA are the ones that have are 'out_of_range'
nr_RR_fail = sum(RR_val == "fail", na.rm = T),
nr_RR_val = sum(RR_val == "pass", na.rm = T),
nr_amp_total = n() - sum(is.na(amp_seq_val)), # here NA are the ones 'not_included'
nr_amp_fail = sum(amp_seq_val == "fail", na.rm = T),
nr_amp_val = sum(amp_seq_val == "pass", na.rm = T),
nr_compound_total = n() - sum(is.na(amp_seq_val)), # here NA are the ones 'not_included
nr_compound_fail = sum(compound_val == "fail", na.rm = T),
nr_compound_val = sum(compound_val == "pass", na.rm = T)) %>%
mutate(perc_qPCR_val = nr_qPCR_val/nr_qPCR_total,
perc_RR_val = nr_RR_val/nr_RR_total,
perc_amp_val = nr_amp_val/nr_amp_total,
perc_compound_val = nr_compound_val/nr_compound_total) %>%
ungroup()## `summarise()` has grouped output by 'tool'. You can override using the
## `.groups` argument.
perc_valperc_val %>% write_tsv('../data/Supplementary_Table_4_precision_values.txt')perc_val %>%
filter(count_group == 'count ≥ 5') %>%
pull(perc_RR_val) %>%
quantile()## 0% 25% 50% 75% 100%
## 0.7402597 0.9491237 0.9625000 0.9808504 1.0000000
perc_val %>%
filter(count_group == 'count ≥ 5') %>%
pull(perc_qPCR_val) %>%
quantile()## 0% 25% 50% 75% 100%
## 0.9000 0.9875 0.9875 1.0000 1.0000
perc_val %>%
filter(count_group == 'count ≥ 5') %>%
pull(perc_amp_val) %>%
quantile()## 0% 25% 50% 75% 100%
## 0.3000000 0.9329323 0.9552239 0.9843137 1.0000000
perc_val %>%
filter(count_group == 'count ≥ 5') %>%
pull(perc_compound_val) %>%
quantile()## 0% 25% 50% 75% 100%
## 0.2714286 0.9047379 0.9310345 0.9527530 0.9830508
first calculate the total nr of validated circ per count group
nr_val = cq %>%
# get set of uniquely validated circRNAs
filter(qPCR_val == 'pass',
RR_val == 'pass',
amp_seq_val == 'pass') %>%
select(circ_id, cell_line, count_group_median) %>% unique() %>%
count(count_group_median) %>%
rename(nr_expected = n)
nr_valthen calculate sensitivity by dividing nr of circ found by total
sens = cq %>%
# get set of uniquely validated circRNAs
filter(qPCR_val == 'pass',
RR_val == 'pass',
amp_seq_val == 'pass') %>%
select(circ_id, cell_line, count_group_median) %>% unique() %>%
# check witch tools have detected these
left_join(all_circ %>%
select(tool, circ_id, cell_line) %>% unique()) %>%
group_by(tool, count_group_median) %>%
summarise(nr_detected = n()) %>%
left_join(nr_val) %>%
mutate(sensitivity = nr_detected/nr_expected) %>%
ungroup()## Joining, by = c("circ_id", "cell_line")
## `summarise()` has grouped output by 'tool'. You can override using the
## `.groups` argument.
## Joining, by = "count_group_median"
senssens %>% write_tsv('../data/Supplementary_Table_5_sensitivity_values.txt')RT-qPCR
cq %>%
select(circ_id, cell_line, qPCR_val) %>% unique() %>%
count(qPCR_val)round(1479/(1479+37), digits = 3)## [1] 0.976
RNase R
cq %>%
select(circ_id, cell_line, RR_val) %>% unique() %>%
count(RR_val)round(1319/(1319+85), digits = 3)## [1] 0.939
round(112/(1319+85+112), digits = 3)## [1] 0.074
amplicon sequencing
cq %>%
select(circ_id, cell_line, amp_seq_val) %>% unique() %>%
count(amp_seq_val)cq %>%
select(circ_id, cell_line, amp_seq_val_detail) %>% unique() %>%
count(amp_seq_val_detail)round(1014/(1014+165), digits = 3)## [1] 0.86
round(337/(1014+165+337), digits = 3)## [1] 0.222
nr of unique circ that pass all val methods
cq %>%
select(circ_id, qPCR_val, RR_val, amp_seq_val, cell_line) %>%
unique() %>%
filter(qPCR_val == 'pass', RR_val == 'pass', amp_seq_val == 'pass')round(957/1516, digits = 3)## [1] 0.631
nr of unique circ that fail all val methods
cq %>%
select(circ_id, qPCR_val, RR_val, amp_seq_val, cell_line) %>%
unique() %>%
filter(qPCR_val == 'fail', RR_val == 'fail', amp_seq_val == 'fail')round(18/1516, digits = 3)## [1] 0.012
nr of unique circ that fail 1 or 2 val methods
cq %>%
select(circ_id, qPCR_val, RR_val, amp_seq_val, cell_line) %>%
filter(!is.na(RR_val), !is.na(amp_seq_val)) %>%
unique() %>%
mutate(val_nr = str_count(paste(qPCR_val, RR_val, amp_seq_val),
pattern = 'pass')) %>%
#filter(val_nr == 0)
#filter(val_nr == 3)
filter(val_nr < 3, val_nr > 0)round(128/1516, digits = 3)## [1] 0.084
nr of cicrc that have at least one NA
cq %>%
select(circ_id, qPCR_val, RR_val, amp_seq_val, cell_line) %>%
unique() %>%
mutate(NA_nr = str_count(paste(qPCR_val, RR_val, amp_seq_val),
pattern = 'NA')) %>%
#filter(NA_nr == 0)
#filter(NA_nr == 3)
filter(NA_nr > 0)round(413/1516, digits = 3)## [1] 0.272
nr of unique circ that have no NAs
cq %>%
filter(!is.na(RR_val), !is.na(amp_seq_val)) %>%
select(circ_id, qPCR_val, RR_val, amp_seq_val, cell_line) %>%
unique()=> recalculate previous numbers based on no NAs all pass
round(957/1103, digits = 3)## [1] 0.868
all fail
round(18/1103, digits = 3)## [1] 0.016
one or two fail
round(128/1103, digits = 3)## [1] 0.116
RT-qPCR
cq %>%
filter(count_group == 'count < 5') %>%
select(circ_id, cell_line, qPCR_val) %>% unique() %>%
count(qPCR_val)round(243/(243+17), digits = 3)## [1] 0.935
RNase R
cq %>%
filter(count_group == 'count < 5') %>%
select(circ_id, cell_line, RR_val) %>% unique() %>%
count(RR_val)round(165/(165+25), digits = 3)## [1] 0.868
amplicon sequencing
cq %>%
filter(count_group == 'count < 5') %>%
select(circ_id, cell_line, amp_seq_val) %>% unique() %>%
count(amp_seq_val)round((63+136)/(63+136+61), digits = 3)## [1] 0.765
round(136/(63+136), digits = 3)## [1] 0.683